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^' ABSTRACT 

<: 

^ [ The Athena MHD code has been extended to integrates the motion of parti- 

I cles coupled with the gas via aerodynamic drag, in order to study the dynamics 

q"'! of gas and sohds in protoplanetary disks and the formation of planetesimals. Our 

nj I particle-gas hybrid scheme is based on a second order predictor-corrector method. 

r^ I Careful treatment of the momentum feedback on the gas guarantees exact con- 

Y^', servation. The hybrid scheme is stable and convergent in most regimes relevant 

2 ! to protoplanetary disks. We describe a semi-implicit integrator generalized from 

c/3 I the leap-frog approach. In the absence of drag force, it preserves the geometric 

i| properties of a particle orbit. We also present a fully-implicit integrator that is 

—^ I unconditionally stable for all regimes of particle-gas coupling. Using our hybrid 

^ I code, we study the numerical convergence of the non-linear saturated state of 

QQ I the streaming instability. We find that gas flow properties are well converged 

QJ^ I with modest grid resolution (128 cells per pressure length rjr for dimensionless 

• I stopping time r^ = 0.1), and equal number of particles and grid cells. On the 

O ' other hand, particle clumping properties converge only at higher resolutions, and 

o: 



% 



finer resolution leads to stronger clumping before convergence is reached. Fi- 
nally, we find that measurement of particle transport properties resulted from 
the streaming instability may be subject to error of about ±20%. 

Subject headings: hydrodynamics — instabilities — methods: numerical — plan- 
etary systems: protoplanetary disks — turbulence 



1. Introduction 

Aerodynamic coupling between gas and solid bodies plays a crucial r ole in the dynam- 



ics of protoplanetary d isks (PPDs) and in planetesimal formation (e.g., ICuzzi et al. 



1993 



Chiang fc Youdinll2009l ). In PPDs, the gaseous disk is partially supported by the radial pres- 
sure gradient, and rotates at sub-Keplerian velocity. On the other hand, the solid particles 
are not affected by the pressure gradient and tend to orbit at Keplerian velocity, resulting 
in relative motion and gas drag. The drag force is characterized by the stopping time tstop, 
where in the absence of other forces, the particle velocity relative to gas would decrease with 
time as exp (— t/tstop)- It is most conveniently parameterized by the dimensionless stopping 
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time, Ts = fitstop, where Q is the Keplerian angular frequency. Tg measures the strength of 
couphng between gas and sohds and depends strongly on particle size and locatio n in the 
PPDs. In the Epstein regime (for partic le size smaller than th e gas mean free path, [Epstein 
19241 ). which is most relevant in PPDs ( IWeidenschiUingl 119771 ) . r^ = psttfl/pgCs, where ps is 
the solid density of the particles, a is the particle radius (assuming spherical shape), pg is gas 
density and c.s is gas sound speed. For the standard minimum mass solar nebular (MMSN) 
model (JHayashil 1 19811 1. r^ = 1 roughly corresponds to meter sized bodies at 1 AU. 



In PPDs, dust grains settle to the disk midplane layer, a nd the interaction betw een gas 
and solids may generate Kelvin-Helmholtz instabilit y (KHI) f lWeidenschilling] Il980l ) and/or 
the streaming instability (JYoudin fc Goodmanll2005l ). Hybrid simulations with both gas and 
solids are necessary to uncover the non-linear dynamics in the dusty midplane layer. In 
addition to the hydrodynamic solver, two ingredients are essential for such hybrid simula- 
tions: First, large number of super-particles are required to mimic the size distribution of 
real solids in the PPDs. Each super-particle represents a swarm (billions or more) of real 
particles with the same phy sical properties. Alth ough a two-fluid approach can be used to 
model the dusty disks (e.g., iGaraud &: Linll2004l ). the particle approach is necessary to ac- 
count for a distribution of solid bodies with different physical properties, and when the solids 
are largely collisionless. The second ingredient is the aerodynamic interaction between gas 
and particles. In particular, the momentum feedback from particles to gas must be included, 
which is essential for the development of KHI and streaming instability. 

Hybrid codes of the kind mentioned above have been developed by several g roups 
jYoudin fc JohansenllJoOTI : I Johansen &: Youdinll2007l : lBalsara et al.ll2009l : iMiniatillJoioh . One 
of the primary goals of this paper is to present the implementation of hybrid particle-gas 
integration scheme into the Athena code, a new grid-based code for compressible magnetohy- 
drodynamics (MHD) based on higher order Godunov methods. A c omprehensive descr iption 



of the implementation and tests of the MHD algorithms are given in lStone et a. 

under lying 

J2009h and 



lydrodynamic solver in Athena is similar to the methods used in 



J2008h. The 



Balsara et al. 



Miniatil (120101) . but quite different from the fin ite difference methods, used by 



Youdin fc JohansenI (120071 ) and I Johansen &: YoudinI (120071 ) (hereafter, YJ07 and JY07 re- 
spectively). Our implementation of the hybrid particle-gas scheme is different from any of 
the previous methods and has several novel features. First of all, in solving the coupled 
equations of gas and particles in hybrid simulations, the problem becomes difficult when the 
particles are strongly coupled to the gas (i.e., r^ ^ 1), so that the coupled equations become 
very stiff. This is especially relevant to submillimeter dust grains. We have developed a 
semi-implicit and a fully-implicit particle integrators that can handle most particle-gas cou- 
pling regimes relevant for PPDs. Secondly, our semi-implicit particle integrator naturally 
generalizes from the leap-frog type integrator, and preserves geometric properties of particle 
orbits. Thirdly, our hybrid scheme conserves linear momentum to machine accuracy. 

We show a suite of test problems to demonstrate the performance of our particle inte- 
grators and the hybrid scheme. In particular, we show that using the fully-implicit integrator 
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is necessary when the particle stop ping time is less than the numerical time step. Further , 
we test and compare our code with JYoudin fc JohansenI (120071 ). I Johansen fc YoudinI ( 120071 ). 
Balsara et al.l (|2009[ ). and iMiniatil ( l2010l ). YJ07 provides two eigen- vectors of the unstable 
modes of the streaming instability. Measuring the linear growth rates of these unstable 
modes and comparing them with theoretical values constitutes a stringent test problem of 
the hybrid code. JY07 studied the non-linear saturation of the streaming instability with 
a set of model parameters. They investigated a variety of physical quantities as measured 
from their simulations. We show that our code reproduces these published results. 



Finally, we present a systematic study of the numerical convergence on the simula- 
tions of the streaming instabilities. Although JY07 performed some experiments to test the 
convergence of their numerical results, they did not carry out a systematic study on this 
subject. We pick two representative runs from JY07, and repeat the simulations with dif- 
ferent grid resolution, and different number of particles in the simulation box. The results 
from our study provides useful insights in understanding the uncertainties of various physical 
quantities as measured from such hybrid simulations. 

This paper is organized as follows. The formalism and the hybrid particle-gas inte- 
gration scheme of our code is presented in §2j In ^ we describe and discuss two implicit 
particle integrators that solves the stiffness problem. We provide a suite of test problems 
to demonstrate our code performance in §H In ^ we systematically study the numerical 
convergence of the streaming instability in the saturation state. We summarize and conclude 
our paper in §Si 



2. Hybrid Particle-Gas Scheme 



2.1. Formalism 



Motivated by the stu dy of the streaming instability in the context of PPDs (e.g., 
Youdin fc Goodmanll2005(). we formulate th e equations with the local shearing sheet ap- 
proximation (JGoldreich fc Lynden-Belllll965[ ). The source terms in this approximation can 
be dropped to study other problems. We choose a local reference frame located at a fiducial 
radius, corotating at the orbital angular velocity Q. The dynamical equations are written 
using the Cartesian coordinate, with x, y, z denoting unit vectors pointing to the radial, az- 
imuthal and vertical direction. In this non-inertial frame, the coupled equations of particles 
and gas read 



dt ^ 



u 



stop 
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+ V ■ {PgU) = 
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In the above equations, pg, Pg denote the mass density and pressure of the gas, q = 
dlnQ/dlnr is the background shear parameter, with q = 3/2 for Keplerian flow, and u, 
V denote velocities of the gas and particles in this reference frame. The subscript "z" in 
equation ([T]) represents the ith. particle. The p article stopping time d ue to gas drag, tgtop, 
depends on particle size and gas flow properties ( IWeidenschillingj 119771 ). Our code is capable 
of dealing with an arbitrary number of different particle species (each particle species has a 
different stopping time), but unless otherwise stated, we assume single particle species with 
constant stopping time throughout this paper for simplicity. In equation 02]), v stands for 
averaged particle velocity in the "fluid element" (weighted by mass), and e denotes the local 
mass density ratio between particle and gas e = Pp/ Pg- This term represents momentum 
feedback from the particles to the gas, written in the form of treating particles as a fluid. 
The particle treatment of feedback term is described in §2.2[ and conservation of total mo- 
mentum is guaranteed. In this paper, we consider non-stratified disks by neglecting vertical 
gravity terms in the equations above (i.e., the Vt^z terms). We also neglect terms associated 
with the m agnetic field in th is paper. They are handled by the underlying MHD integrators 
in Athena ( IStone et al.ll2008[) . An isothermal equation of state for the gas is used throughout 
this paper, with P 



Pgcl 



f 



Our goal is to perform the local shearing box simulations (JHawley et al.lll995l ). where 
the radial boundary condition is periodic with additional shear to account for differential 
rotation. Therefore, it is not appropriate to include radial pressure gradient directly, which 
is inconsistent with the periodic boundary conditions. Alternatively, one can replace the 
pressure gradient by a constant radial force acting on the gas F = 2'qvK^x, pointing outward. 
The quantity rjvx measures the amount by which the gas (azimuthal) velocity is reduced from 
the Keplerian value due to the radial pressure gradient. In our code, instead, we find it more 
convenient not to modify the hydrodynamic integrator, but to add a constant radial force 
on the particles, pointing inward. Our treatment is mathematically equivalent to the effect 
of a radial pressure gradient, but both particle and gas (azimuthal) velocities are shifted to 
slightly larger value, by rjvx. 



An orbital advection scheme (JMassetll200d : I Johnson et al.ll2008l ) has been implemented 
in Athena, which takes the advantage of the fact that the above equations can be split into 
two systems, one of which corresponds to linear advection operat or with background flow 
velocity {qVtx)y and the other involves only velocity fluctuations ( IStone fc Gardinerll2010l ). 
We use the orbital advection scheme in all our simulations, not only because it is it is faster, 
but also more accuratcl- The same technique can be implemented to the particles, with 



^ The drag between gas and solids dissipates energy and generates heat. It is ignored in the isothermal 



gas. 



^This scheme is used only in 3D shearing box simulations (e.g.. lBai fc StonduOlOa ). In this paper, where 
we perform 2D simulations in the radial-vertical plane, the orbital advection scheme is not used. 
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equation ([T| 



replaced by 

dt 
{qfix)y, and u' 



2« - VVk)^x - (2 - q)ii^ny - n^zz 



f • 



u 



t 



(4) 



stop 



where v' = v — [qVtx)y, and u' = u — {qVtx)y. An particle advection step of is then carried 
out in parallel with the orbital advection of gas. Also note that we've included the effect of 
gas radial pressure gradient in the first term on the right hand side of equation (jll). 



2.2. Predictor-Corrector Scheme 



The MHD integrator in Athena adopts a directionally unsplit, higher-order Godunov 
method, which conserves mass, momentum and energy (when applicable) to machine accu- 
racy. Our goal is to develop a particle-gas hybrid integrator that is also conservative and at 
least second order accurate. Two MHD integrators have been ir nplemented in the Athena 
code, including the corner ed transport upwind (CT U) integrator (IStone et al.ll2008l ) and the 
van-Leer (VL) integrator (jStone fc Gardinerll2009l ). Our parti cle scheme is combined with 
the CTU integrator, which is more accurate and less diffusive (jStone fc Gardinerll2009l ). In 
our implementation, the mornentum (and energy) feedback to the gas is treated as source 
terms (different from lMiniatill2010l ). while the evaluation of transverse flux gradient in the 
hydrodynamic solver is unchanged. The gas continuity equation is automatically handled 
in the Godunov scheme, with no modifications to the code needed. The hybrid particle-gas 
scheme adopt a predictor-corrector approach, which is described below. 

To demonstrate the numerical algorithm, we rewrite the coupled particle-gas momentum 
equations. For the gas momentum equation, we simplify the left hand side of equation ([3]) 
to a Lagrangian derivative, since we do not modify the calculation of the flux gradients in 
the CTU integrator. We use /(aj, v) and g{x^ u) to denote source terms for the acceleration 
of particles and gas respectively, due to forces other than the drag force (e.g., Coriolis force, 
tidal force and the global pressure gradient). The formalism in §2. II can be summarized as 

v.; — U 
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We do not distinguish between (it, v) and {u', v') here because it does not affect our de- 
scription of the numerical algorithms. Our predictor-corrector scheme, which integrates the 
coupled equations from time step t^^'^ to t^"'~^^\ can be illustrated as 
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In the above, ( !7a|) and (17cl) are generalizations of the predictor and corrector step of the 
MHD integrator, and particle feedback to gas is expressed as Ap^'^" andAp^™'' for the two 
steps respectively. Their expressions are given in the following paragraphs. W is the weight 
function for interpolation (see next paragraph). fl7bp represents the particle integrator, which 
we will discuss in detail in ^ Note that in the bracket on the right hand side of this equation, 
particle quantities Xi and Vi can be combinations of step (n) and step (n + 1) quantities, 
depending on the particle integrator. 

The calculation of the drag force experienced by particles requires interpolation of grid 
quantities to the particle location on the grid. The interpolation scheme is described by 
the weight function W{x — Xi). For consistency, the same interpolation scheme is used to 
distribute the feedback from individual particles to the gas grid points, as in equation ([7]). 
In order avoid spurious numerical artifacts in the hybrid scheme, the weight function W 
should satisfy certain constrains. It should be continuous over the computational domain 
to avoid sharp transitions , and W{Ax) should be non- negative for any Ax to reduce noise 



( jYoudin fc Johansenll2007l ). The interpolation should be accurate enough to minimize errors. 
In particular, interpolation error should not be much worse than the error from the spacial 
reconstruction of the MHD integrator (we use the third order piecewise parabolic method). 
Finally, interpolation is time consuming, so the s cheme should be as efficien t as possible. We 



have compared three interpolation schemes (cf lBirdsall fc Langdoru l2005l YJ07), namely, 
cloud-in-a-cell (CIC), triangular-shaped cloud (TSC) and quadratic polynomial (QP). The 
CIC scheme is simple but inaccurate and noisy, the QP scheme is the most accurate but 
not continuous. Similar to JY07, we choose the TSC interpolation scheme throughout this 
paper. 

In the predictor step, the momentum feedback from individual particles is calculated 
from direct force evaluation, multiplied by half a time step. However, when the particle 
stopping time is small, the drag force diverges, and the error in the velocity calculation and 
interpolation can be substantially amplified. We note that in the absence of other forces, 
the particle velocity would approach gas velocity as exp (— t/tstop)- Therefore, we modify the 
predictor step momentum feedback (from individual particle "z") into 

^ prcd ^ m,{vi-u) h 
' max (tstop, h) 2 

where mi is the particle mass. We do not use the exponential expression because it is derived 
by assuming drag force only. Our treatment is the same order accurate as the exponential 
expression (first order) when tstop < h, and when tgtop > h, it ensures exact force balance at 
equilibrium statqj. The individual particle feedback is then distributed to neighboring grid 
cells. As a source term, the momentum feedback is divided by gas density and is added to 



•^For example, when testing the hnear growth rate of the streaming instabihty (in Q, one starts from 
the NSH cquihbrium. Our treatment aUows NSH equihbrium to be satisfied exactly. 
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the left and right states of priin itive variables (gas velocity) between step 1 and step 2 of the 



CTU algorithm as described in IStone et al.l ( 120081 ) 



For the feedback calculation in the corrector step, note that the particles have already 
evolved for a full time step, we can calculate the momentum difference of individual particles 
between the two steps (n) and (n + 1). Let Fc denote forces experienced by particles other 
than the drag force. Since Fc is generally non-stiff (e.g., Coriolis force), we can obtain 
momentum feedback from particle i to be 

where F^^ is evaluated at x^'^~^^^'^\ the midpoint between x^"'^ and x^"'~^^\ which ensures 
second order accuracy. This treatment is conservative and guarantees exact momentum 
conservation. Moreover, it avoids the potential stiffness due to the direct evaluation of the 
drag force. To distribute the feedback to grid points, we again take the force location at 
^(n+i/2)^ as indicated in equation (17c|) . The corrector step feedback is added to the end of 
the gas integrator. 

If one considers the disk thermodynamics (not treated in this paper), the heat generated 
by friction needs to be deposited to the energy equation of the gas, with energy dissipation 
rate (per unit volume) 

£ = pp{v-uY/tstop ■ (10) 

The heating term is contributed from the work done by the pressure gradient and is typically 
only a small amount compared with the disk thermal energy budget (see YJ07 for more 
details). Numerically, potential stiffness problem also exists since tstop is present in the 
denominator. In practice we rewrite the energy deposition rate from particle i to be Si = 
{Api)Hstop/'m'ih'^ , where m is particle mass and Api is taken from equations (|8]) and ([9]) for 
predictor and corrector steps. 

The overall accuracy of our hybrid scheme is second order, less than the Pencil code, 
which is a finite-difference code with higher order accuracy in smooth flows. However, our 
code is fully conservative, and it does not need to be stablized by artificial hyper-viscosity. 
It also all ows the development of implicit particle integrator (see next section) much easier. 



Recently, iBalsara et al.l (J2009[ ) has proposed a similar predictor-corrector hybrid scheme for 
particle-gas dynamics, but approximation is made in the corrector step (see their equation 
(16) and the d iscussion that follows) and is less accurate when dust dominates local density. 



Miniatil ( 120101 ) has described another hybrid scheme that is fully implicit, but assumes only 
one particle species. Our predictor-corrector scheme is simpler than these other approaches, 
allows an arbitrary number of particle species, while we will show in ^and §3] that our code 
performance is at least comparable to all these codes. 

One disadvantage of our predictor-corrector hybrid scheme is that it is intrinsically 
explicit, and may cause numerical instability when the coupled equations are stiff. Here, the 
stiffness is caused by the parameter e, the local particle to gas mass ratio. The design of the 



particle integrators assumes that particles move in the gas velocity field. In the regime where 
4top ^ h, our semi-implicit and fully-implicit schemes force particles to be strictly coupled 
to the gas. However, when e ^ 1, gas is expected to follow the particles. This situation may 
result in unphysically large growth of particle and gas velocities, making the hybrid code 
unstable. More specifically, we define the stiffness parameter 

X = ^ Cfc/i/ max (tstop,fc, h) . (11) 

k 

In the above expression we have generalized our analysis to include a size distribution of 
particles, and subscript "/c" labels different particle types. Our hybrid scheme can become 
unstable when x exceeds order unity. Our experiments show that the threshold value of x 
is about 3 — 5 (depending on the actual problem). 

One way to remedy the stiff particle mass loading problem discussed above is to make 
the overall hybrid scheme implicito. This is relatively easy to do if all particles h ave a single 



size, and the treatment will be similar to a two-fluid approach (JMiniatill20ld ). However, 
the two-fluid treatment is not easily generalized to a size distribution of particles. This is 
because particles with different sizes are (indirectly) coupled to each other via interactions 
with the gas, and evolving the coupled equations implicitly requires solving an inverse matrix 
of rank A^ -|- 1 at each grid cell, where N is the number of particle size bins. The matrix 
is non-diagonal due to the coupling U- As A^ gets larger and larger, the algorithm becomes 
more and more complicated and computationally prohibitive. 

Alternatively, in regions where x exceeds certain thresholds, one can effectively increase 
the value of tgtop for all local particles to bring down the value of x, which guarantees stability, 
as is adopted in the Pencil Code (A. Johansen, private communication, 2010). Physically, 
this means in dense particle clumps, particles can move more freely and are less affected by 
gas drag, while just sufficient momentum feedback is added to the gas for it to follow the 
motion of particles without causing any numerical instability. We have implemented this 
technique that enforces x < 3 in our code, although in practice, this feature is turned off 
due to the reasons below. 

Fortunately, in the context of PPDs, the overall dust to gas mass ratio is about one 
percent. Therefore, in an average sense, x is much less than unity, x can become larger 
when dust grains settle towards the disk midplane. The typical time step in the simulations 
is about Qh = 10~^ — 10~^. For dust grains with stopping time Ts < 10"^, one can easily 
show that even very weak turbulence is able to make these solids suspended in the disks, 
keeping x relatively small at disk midplanqj. Concentration of particles can raise x locally. 



** Another approach is to artificially reduce the momentum feedback by max [x, 1] in the predictor step 
(not in the corrector step, otherwise momentum conservation would break down). This approach reduces 
the accuracy of the code to first order, but improves the stability. 

^An example of such a matrix is given in the Appendix A of iBai fc Stond ( 2010a ). 

^With the standard a prescription ([Shakura fc Sunvaevlll973[ ). the average value of x in disk midplane 
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In the case of the streaming instabihty, however, particle concentration is most efficient for 
marginally coupled solids r^ ~ 1. Clumping of these particles only raises x very slowly 
according to equation fITT]) . In our simulations of the SI in §3], the maximum value of x 
never exceed the threshold x = 3 (see §5.21 for more details). We emphasize that the value 
of X only determines stability, but does not constrain the maximum particle density. In 
some of our stratifi ed disk simulations with large solid abundance (5-7 times super-solar. 



Bai fc Stonell2010bl ). we do observe the numerical instability and have to reduce the time 



step in our calculations. 



3. Particle Integrator 

In this section, we describe in detail the particle integrator that have been implemented 
in our particle-gas hybrid scheme [i.e., equation (]7bp ]. The overall problem for the particle 
integrator is 

^ = ^^, ^ = a[^,x,n('^+V2)(^)]^ (12) 

where a denotes the acceleration due to all the forces, including the gas drag, and following 
the convention of equation ((Tj). As suggested in equation (l7bl) . we use half time step gas 
velocities n'^""'"^/^) to avoid tracking the evolution of gas, which also ensures second order 
accuracy. For the sake of simplicity, in the remaining of this section, we will drop the 
superscript ("+1/^) in the gas quantities. 

As we have noted before, the gas drag term becomes stiff for strongly coupled particles. 
We have developed two second-order implicit particle integrators. Each integrator has its 
own pros and cons. In our code, we allow different species of particles to be pushed by 
different integrators, which enable us to integrate particles with any stopping time while 
maintaining the geometric properties of particle orbits. 



3.1. A Semi-implicit Integrator 

Our first integrator is a semi-implicit scheme based on the Crank-Nicholson method. 
The basic algorithm is the following. 



X =x 



(n) 



/iv("V2 



-y(" + l) 



v(") + /ia[(t;(") + 'j;("+^))/2, x' , u{x')] , (13) 



can be estimated by % '^ Z(niax [flh, Ts]/a)^''^, where Z is the ratio of dust surface density over gas surface 
density, typicaUy Z ~ 0.01. Thus fairly weak turbulence a — 10~^ is sufficient to keep x well below unity. 
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where x' = x^"'^ + hv^^'' /2 is the predicted particle position at t("+^/2), and is used to evaluate 
the stopping time and the gas velocity. This scheme is semi-implicit in the sense that the 
velocity update depends on velocities in both step [n) and step (n + 1). Converting to explicit 
form, we find 

^(n+i) _ „(n) = hA-'aiv^^'lx', u{x')) , A = 1 - -^ , (14) 



2dv 



where the Jacobian da/dv is evaluated at x', with 



da 
dv 





(15) 



We always use the orbital advection algorithm ( IStone &: Gardinerll201Cll ). therefore the (2 — g) 



factor in the matrix element. After calling the particle integrator, we need to shift the 
particle positions by an amount —qVtxhy. In practice, we replace x by [x*^"-* +x*^"'"^^'']/2. We 
emphasize that using the particle advection scheme is important for improving the accuracy 
of the particle integrator, and especially, preserve the geometric properties of particle orbits 
(see below). 

This semi-implicit integrator has close analogy to the Boris integrator due to the sim- 
ilarity between Coriolis force and Lorentz force. It is essentially a Leap-Frog integrator, in 
the form of "Drift-Kick-Drift" (DKD). In ^U we will see that in t he limit tstnp = oo, this 



integrator preserves geometric orbital properties exactly. Recently, iQuinn et al.l ( l2009l ) has 
proposed a symplectic particle integrator for Hill's equations, which is essentially in the form 
of "Kick-Drift-Kick". In their work, particle advection was not implemented, making their 
formula slightly more complicated than ours. 

This integrator, although is implicit, can still cause problems when the particle stopping 
time approaches zero (i.e., when tstop '^ h). In this limit, the position update reduces to 
^.{n+i) _ ^(n) j^fiiii^x')^ which is indeed second order accurate. However, the velocity update 
reduces to i?^"^ + i;('"+^) = 2u{x'). If the initial velocity difference between particle and gas 
is large, then particle velocity will oscillate around the gas velocity without damping. More 
seriously, the evolution of gas velocity may even amplify the velocity difference between 
particle and gas, causing a runaway. Our experiences indicate that this integrator works 
safely at least for t > 0.2/i. One can construct other similar second order semi- implicit, 
Crank-Nicholson type schemes. Nevertheless, the only other possibility one can achieve in 
the tstop — ^ limit is t;('") + t>("-+i) = ■u(x'''"^) + it(a;("+-^^), which is prone to the same problem. 
Also, other schemes no longer maintain geometric properties of particle orbits. 



3.2. A Fully-implicit Integrator 

The apparent weakness of any types of the semi-implicit method leads us to develop 
an absolutely stable, second order integration scheme. To be absolutely stable, we demand 
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the method to be fully implicit, that is, the velocity update to step (ra + 1) depends only on 
velocity at step (n + 1). For a simple ordinary differential equation dy/dt = f{y), a second 
order fully-implicit scheme can be constructed as 



yin+l) = yin) + /,/(^("+l))/2 + /^/[y("+^) - /l/(l/("+^))]/2 . 



(16) 



Other second order fully-implicit schemes are possible, but we will stick to this specific form. 
The implementation of this method to the particle integrator is illustrated as follows 
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where we have omitted the fluid argument ii('"+i/2) for conciseness, and we assume that u is 
evaluated at the same position {x' or a;^'"^) as in a. It can be easily shown that in the limit 
^stop = 0, the above equations demand that i)("+^) = u{x'). Although this is not second 
order accurate, the position update is indeed second order. Therefore, we achieve second 
order accuracy with absolute stability. 

The main work is to update the velocity, and to the second order, we find 
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where the subscript "0" means to evaluate the Jacobian at a;'^"^ "1" means to evaluate at 
x' (note that the stopping time can depend on position in the most general case). The 
inverse matrix A~^ can be evaluated analytically without any trouble, since it involves only 
the inversion of a 2 x 2 matrix. Alternatively, it suffices to evaluate the inverse matrix up 
to first order in h, except for terms containing 1/tstop, where all orders of h should be kept. 
The result is 
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From equation (fT9|) . we see that as t^top — > 0, A~^ — )> b~^ ~ 211^^^/11^ — ?> 0. Therefore, 
this integrator is stable for any tstop > 0. One disadvantage of this integrator is that it 
no longer preserves geometric properties of particle orbits as tstop — ;• C)0. In fact, like the 
explicit method, any fully-implicit method always fail to preserve such geometric properties 
(see §4.ip . Therefore, we will in most cases use the semi- implicit particle integrator, while 
we switch to this fully-implicit integrator for particles with tstop ^ h. 
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An alternative way of integrating th e strongly coupled particle s is to use the "short fric- 



tion time" approximation introduced by iJohansen fc Klahrl (120051 ) . In this approximation 



particles can be considered as being carried with the gas, while maintaining a small drift 
velocity due to some external force Fp. Meanwhile, the gas also feels some external force Fg 
other than the drag. Fg is in general different from Fp, since particles don't have pressure 
support, and generally don't feel the magnetic field. In this approximation, particle velocity 
is assumed to be always equal to the termination velocity 



''term I 



u{x)+U,oAFp-Fg) . (20) 



One can then easily integrate particle orbit based on this expression to any order of accuracy. 
We note that leading truncation error of this approximation is of the order max[(tstop| Vm|)^, 
^stopl^-^Mstop^l^'^l]' "where F = Fp — Fg. Roughly speaking, the short friction time 
approximation is applicable for stopping time tstop ^ min(r2~^, iVul"^). Converting the 
above into a second order integrator further introduces truncation error of the order {Qhy. 
Therefore, this integrator may perform equally well as our fuUy-i r nplicit integrator only 



when tstop < h. The original implementation of iJohansen fc Klahrl ( 120051 ) did not include 
the momentum feedback, thus can not be used to study the SI. However, it can be extended 
to include feedback as described in §2.21 Nonetheless, we do not implement this integrator 
in our code since we have the fully-implicit integrator at hand. 



4. Code Tests 

4.1. Epicycle Test 

We begin by examining the performance of the particle integrators in the weak coupling 
limit, i.e., tstop = oo, so gas drag does not enter the problem. From the particle equation of 
motion ([1]), the particle trajectory follows an epicyclic orbit: 

xit) = Acosut , y(t) = Asinut , u = a/2(2 - q)Q . (21) 

u 



where x{t),y{t) denote radial and azimuthal directions relative to domain center, and A is 
the radial amplitude of the oscillation. Epicyclic motion conserves total energy 

E = hx^ + f) - qn^x^ = (2 - q)n^A^ . (22) 

We integrate a test particle and follow its epicyclic orbit using different particle in- 
tegrators and examine the particle trajectory and energy conservation. In particular, we 
consider three 2nd order integrators, namely, the semi-implicit and fully-implicit integrators 
introduced in ^ and we add an explicit integrator based on the modified Euler method 
for comparison purpose. In Figure (H we see that the semi-implicit method conserves total 
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energy exactl}U. The particle orbit is closed, and the truncation error exhibits as a phase 
shift relative to the analytical solution. The phase error diminishes as h"^ (not shown in the 
figure), as expected. Note that in order to make the errors significant, we have chosen a 
rather large time step h = 0.4^2"^. In typical numerical simulations, the time step and thus 
phase error is much smaller. From Figure [H we also see that the test particle gains energy 
and drifts away from the guiding center when using the explicit method, while it loses energy 
and gyrates into the guiding center when using the fully-implicit method. This result is quite 
general. 

This test demonstrates that the semi-implicit method preserves the geometric properties 
of particle orbits, and moreover, it is efficient because it evaluates the drag force (therefore 
interpolation) only once per step. Therefore, in most applications we prefer to use this 
integrator. 



4.2. Particle-Gas Deceleration Test 

In the second problem, we consider the motion of a collection of particles in a uniform 
gas. The spatial distribution of the particles is uniform. Gas drag and feedback are included 
in the test. We work in the frame where the total momentum in particles and gas is zero. 
Let Wq be the particle initial velocity, then all the particles evolve as 

where e is the overall particle to gas mass ratio, and S is the distance a particle has traveled. 

In our first test, we choose e = l,tstop = 2,u7o = 1, and gas with density Pg = 1. 
We evolve the system to tg = 1 with constant time step h using the semi-implicit and 
the fully-implicit integrators. In Figure [2^ we plot the position error AS" at tg = 1 as a 
function of h. It is clear that the cumulative position error scales with h"^, indicating that 
our particle-gas hybrid code has achieved second order accuracy. Figure [2]3 shows the total 
momentum density at tg = 1 as a function of h. We see that our hybrid scheme conserves 
linear momentum within machine accuracy (to the level of 10~^^). The total momentum 
at te = 1 increases with decreasing h, which refiects the accumulation of round-off error 
with increasing number of time steps. For uncorrelated round-off errors, one would expect 
/i~^'^ scaling. The slightly steeper slope reflects a certain level of correlation in the round-off 
errors. 

Next we test the integrators in the stiff regime. To do so, we flx the time step h = 1, 
and choose tstop to be less than h. The other parameters are the same as in the previous test 



^This result is valid only when the orbital advection algorithm is used. Otherwise, similar to the leap-frog 
method, the numerical Hamiltonian is time-dependent and energy is not conserved, but oscillates around 
the true value (if the time step is not constant, the energy may even gradually deviate from the true value). 
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Fig. 1. — Results of the epicycle test, using f] = 1, g = 3/2, and one test particle with 
initial amplitude A = 0.4, at fixed time step h = 0.4. Left: Total energy of the test 
particle integrated using explicit (dashed), semi- implicit (solid) and fully-implicit (dash- 
dotted) methods. Right: Particle radial orbit integrated using the semi- implicit method 
(solid) compared with analytical solution (dashed). 
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Fig. 2. — Error in position (a) and total momentum density (b) as a function of time step 
length h in the particle-gas deceleration test. We choose the particle to gas mass ratio to 
be 1, and the initial velocities of particles and gas are both 1, but in opposite directions. 
The particle stopping time is set to tgtop = 2. Errors are measured at t = 1. Red squares: 
semi-implicit integrator; Black circles: fully-implicit integrator. 
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Fig. 3. — Numerical evolution of particle velocity in the stiff regime of the particle-gas 
deceleration test. The time step is fixed to h = 1. Panels (a) and (b) show the results from 
the semi-implicit integrator and fully-implicit integrator respectively. Black circles connected 
by red solid lines: tstop = 0.2. Black squares connected by blue solid lines: tgtop = 0.02. In 
panel (b), we further show expected evolution of velocity in the two cases: Red dashed line 
for tstop = 0.2 and blue dash-dotted line for tstop = 0.02. Note that in panel (a) we use a 
linear scale while in panel (b) we use a log scale. 
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(e = 1, t^o = !)• We evolve the system to tg = 10, using tgtop = 0.2 and tgtop = 0.02. In Figure 
Owe show the time evolution of particle velocity from the semi-implicit and the fully-implicit 
integrators. We see that with the fully-implicit integrator, the particle velocity (thus the gas 
velocity, due to momentum conservation) rapidly drops to zero, and smaller tgtop leads to 
faster damping [actually v^'^'^^^ ~ 2(tstop/^)^'i^^"'']- Because tstop is not resolved, the damping 
rate is slower than theoretical values, but still is rapid as relative velocity v drops several 
orders of magnitude in each time step. With the semi-implicit integrator, however, the 
particle velocity undergoes damped oscillation. The smaller tgtop is, the slower the particle 
velocity damps [actually t)("+^) ~ —(1 — 4tstop/^)^ ], opposite to the trend found in the 
fully- implicit integrator. In the limit tgtop = 0, the particle velocity never dies away. The 
behavior of the semi-implicit integrator in this regime exactly resembles a damped harmonic 
oscillator, which makes it difficult for particles to get rid of extra relative velocity with respect 
to gas. This situation is problematic because if the gas itself undergoes some other oscillation 
with a similar frequency, the amplitude of the oscillator may be quickly amplified rather than 
damped due to (numerical) resonant interactions. Therefore, for stability considerations, we 
use the fully-implicit integrator when tgtop ^ h. 



4.3. Streaming Instability in the Linear Regime 

The most stringent test of our code is the streaming instability linear growth rate test 
using eigenmodes provided by YJ07. These eigenrn odes are built from the Nakagawa-Sekiya- 



Hayashi (NSH) equilibrium Nakagawa et al.lll986t also see equation (7) of YJ07]. The NSH 



equilibrium refers to the equilibrium between solids and gas in unstratified Keplerian disks, 
in which gas is partially supported by a radial pressure gradient. Establishing the NSH 
equilibrium numerically requires exact balance among all force terms. With the help of 
the particle advection scheme as well as careful treatment of predictor and corrector step 
momentum feedback (see [|2] and ^ for details), we are able to establish the exact NSH 
equilibrium in our code (i.e., net force on both the gas and particles is zero to machine 
precision). The eigenmodes of the streaming instability are characterized by dimensionless 
wave numbers K^ = kxrjVk/^ and Kz = kzrjVk/^- Taking r/ffc = O.OScg and setting the box 
size to be one wavelength, we have K^ = 0.05(cs/f2)(27r/La;), and similarly in the z direction. 
We fix the box size to be L^ = L^ = 2, therefore, K^ (and K^) sets the sound speed. We 
construct eigenvectors of density and velocity perturbations using formula (10) and Table 1 
in YJ07. We use one particle per cell and each particle is located at cell center initially. In 
order to generate particle density profile of p^ = 1 + A cos k^x cos kzZ (we take A = 10~^), 
we shift particle positions in the radial (x) direction, with the amount of shift proportional 
to cos kzZ. 

The test suite in YJ07 aims at measuring the numerical growth rate of the streaming 
instability eigenmodes as a function of grid resolution. It consists of two problems, both using 
Ts = 0.1 particles. The test "linA" has particle to gas mass ratio e = 3 and K^ = Kz = 30. 
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Fig. 4. — Measure d growth rate of the s t reami ng instabihty eigenmodes "hnA" (a) and "hnB" 
(b) adopted from lYoudin fc Johanseru (120071 ) using the semi-imphcit integrator. Sohd hne 
marks the theoretical growth rate, and dashed lines with symbols label the measured growth 
rate as a function of grid resolution (number of cells per wavelength). Asterisks correspond 
to evolution using fixed CFL number 0.8, diamonds correspond to fixed time step (set by 
CFL=0.8 with 128 cells). 
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The predicted growth rate is s = 0.4190204^2. The second test "hnB" has e = 0.2 and 
Kx = Kz = Q. This test is more challenging because the mode grows very slowly, with 
s = 0. 015476^1^ We consider grid resolutions of 16, 32, 64 and 128 cells per wave length. 
Figure m shows our test results using the semi-implicit integrator. In each test, we conduct 
two set of runs, one with fixed Courant-Friedrichs-Lewy number CFL=0.8, and the other 
with a fixed time step (or CFL=0.1, 0.2, 0.4, 0.8 for each resolution). For reference, the time 
step with lowest resolution (16 cells per wavelength) and CFL=0.8 is Qh = 2.6 x 10~^ for 
"linA" and Qh = 1.3 x 10^'^ for "linB". Both semi-implicit and fully-implicit methods are 
considered, and the results are almost identical (since the time step is much smaller than 
^stop)- We see that our code converges very well to the predicted growth rate in the run 
"linA". Small time step helps with convergence. For the more challenging 'linB" test, we 
get slower convergence. In particular, the mode does not grow for 32-cell resolution using 
CFL=0.8. Better temporal resolution again improves the performance. 

These test results show that with CFL=0.8, about 64 cells are needed to see growth 
(e.g., test linB), although 16 cells seem to be sufficient to capture the most unstable modes 
(e.g., test linA). Viewed from the results, our measured linear growth rates are si milar to 



(and even closer to semi-analytic values than) those of iBalsara et al.l ( l2009l ) and iMiniati 



(120101 ). who employ similar MHD code but different hybrid schemes. Our method is less 
accurate than that presented in YJ07, which is not surprising given the fact that our code 
is second order accurate (both spatial and temporal) while the Pencil code used by YJ07 is 
sixth order spatial and third order temporal accurate for smooth flows, such as in this test. 

The test problems "linA" and "linB" provided in YJ07 both adopt relatively large 
particle stopping time Tg = 0.1. Since we will perform SI simulations on a size distribution 
of particles with stopping time down to 10"^, additional code test is needed to justify the 
ability of our code. In Appendix [Aj we provide two more linear growth tests for particles 
with Tg = 10~^ and Tg = 10~^. Numerical convergence to the analytical growth rate is again 
achieved, although slightly higher resolution is required to reach the convergence. 



5. Convergence of Streaming Instability 

As an application of our hybrid code, we study the streaming instability in the non-linear 
regime. Such calculations have been performed and comprehensively analyzed by JY07. In 
this paper, we focus on the numerical convergence of physical properties in the saturated 
state, which was not fully explored in JY07. Understanding the convergence properties of 
hydrodynamic simulations of the streaming instabilities is important before adding more 
complicated physics (such as MHD), and interpreting the reliability of such simulations. 

We study the numerical convergence using 2D simulations in the x — z (radial-vertical) 



*Also, due to a smaller value of K, the sound speed is smaller than the "linA" test, thus larger time step. 
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plane in which both the grid resolution and the number of (super-) particles vary. By using 
2D simulations, we can explore a wide range of resolutions that is not possible in fully 3D 
simulations. Moreover, the results of 2D and 3D simulations are similar (JY07), thus one 
expects the numerical convergence properties of 2D simulations are similar to those in 3D 
as well. 

The two most significant effects of the streaming instability are the concentration of 
particles into dense clumps, and the generation of turbulence and/or waves. The former 
can be characterized by the probability distribution function (PDF) of particle densities, 
while the latter can be characterized by the turbulent velocity, diffusion coefficient and 
momentum flux. We investigate the numerical convergence of these two aspects in §5.21 and 
§5.31 respectively. 



5.1. Simulation Setup 

We perform simulations of the streaming instability using an isothermal equation of 
state and neglecting vertical gravity. In this case, the properties of the streaming instability 
is largely determined by two parameters, the particle stopping time Tg and the overall mass 
ratio between particles and gas e. The drift velocity rivx does not serve as an independent 
variable, but it sets the length scale of the problem, that is, all length scales can be measured 
in units of rjr = tivk/^- The isothermal sound speed Cg is also a mode l parameter, but it is 



much less important because the gas motion is almost incompressible (JYoudin fc Goodman 



20051 . JY07). We adopt rivx/cs = 0.05 throughout this paper, which roughly matches the 



value expected at lAU in the MMSN model. 

As shown by JY07, there are two basic "modes" to the non-linear saturation of the 
streaming instability. For marginally coupled and larger particles (r^ > 1), the instability 
develops for a wide range of e, around 1. In the saturated state, particles concentrate into 
long stripes, which are mostly aligned with the vertical direction and are slightly tilted in the 
radial direction. For tightly coupled particles (r^ <C 1), prominent instability develops only 
when e > 1, and the turbulent state consists of large voids with narrow particle streams. 
The turbulence is much weaker than the previous case (see Figure 2 and Figure 5 of JY07 for 
each of the two "modes", also see our Figure [5] below) . We have also performed a series of 
simulations with a wide range of parameters (r^ and e), and have reproduced all the results 
in JY07. We use the semi-implicit integrator in all the tests in this section, and results from 
the other two integrators are similar, since the time step is at least 10 times smaller than 
tstop in all test runs. 

To study convergence, we choose two test problems that are representative of the two 
saturation "modes" described by JY07 Q The flrst problem is the same as their run AB. 



'We have also studied the non-Unear behavior of the SI for more strongly coupled particles with different 
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The run parameters include: Tg = 0.1, e = 1.0, L^ = L^ = 2rir. The second problem uses 
parameters from their run BA, with r^ = 1.0, e = 0.2, L^, = L^ = 40?7r. We choose the 
"standard" grid resolution to be 256^, and A^pc = 9 particles per grid cell. The typical time 
step is given by /i ~ CF'L[Lx/N)/cs where N is the number of grid cells in one direction. 
For the standard grid resolution and CFL = 0.8, we have Vth ?« 3 x 10~^ for run AB, and 
Vth pa 6 X 10~^ for run BA. We then conduct a series of runs for each problem. First, we 
change the grid resolution from 64^ to 2048^, with the number of particles per cell fixed at 9. 
Note that increasing spatial resolution at fixed CFL number is accompanied by the increase 
of temporal resolution, which is also important for driving to convergence, as we discussed 
in §4.31 Secondly, we fix the grid resolution to be 256^, and change the number of particles 
per cell, A^pc = 1, 4, 9, 16 and 25. In each of these runs, particles are initialized with random 
positions within the simulation box, with velocities taken from the NSH equilibrium. All 
these runs are tabulated in the first three columns of Table [H where each variation of the 
AB or BA runs is assigned with a run number. Our fiducial runs are AB3 and BA3, while 
our runs AB9 and BA9 have exactly the same numerical parameters as that in JY07. 

The AB runs are fully saturated after about 16^2^^ and the BA runs do not saturate 
until about 160^2"^. We choose the saturation time for the two runs to be some later time 
Ts = 30r2~^, Ts = 300r2~^, after which we start collecting data from the simulations 
and perform measurements. In Figure O we show images of the particle density with different 
grid resolution (64^, 256^, 1024^ from top to bottom) for runs AB (left) and BA (right) at 
saturation (t = T^). Our standard runs (middle panels) resemble the last panels of Figure 4 
and Figure 5 in JY07. For the BA run series, we see that the bulk patterns of particle density 
do not vary much with grid resolution, although more and more detailed features are revealed 
in the higher resolution simulations. For the AB run series, all three grid resolutions show 
cavitation and particle streams as described in JY07, however, the scale of particle clumping 
becomes smaller as one uses finer resolutions. A resolution of 256^ appears to be necessary 
to capture the prominent features in the particle density pattern. The results of these runs 
will serve as the starting point of our more quantitative study of numerical convergence in 
the next two subsections, using physical quantities averaged over the period between Tg and 
T„ where t]^^^ = 80^"^ for AB runs and T^^^^ = 1500^"^ for BA runs. 



5.2. Particle Concentration 

Probably the most important property of the streaming instability is its ability to con- 
centrate particles. In JY07, it was shown that the maximum particle density resulting from 



resolutions. Briefly speaking, finer resolution is needed to capture the SI as Ts gets smaller. Particle clumping 
becomes weaker as t, decreases and reduces to modest overdensities as Ts < 0.01 (with e > 1). The results 



are in qualitative agreement with the linear analysis bv lYoudin fc GoodmanI (|2005[ ). In this paper, however, 
we focus on the non-linear SI tests in JY07 for conciseness, the main conclusions are also applicable to the 
case with more strongly coupled particles. 
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Fig. 5. — Snapshots of particle densities in the x — z plane of run AB (left) and BA (right) 
at saturation of the streaming instability T = Tg with different grid resolutions: 64^ (top), 
256^ (middle) and 1024^ (bottom). All simulations use 9 particles per cell. Particle densities 
are shown with log scale, ranging between 0.1 to 10 of the gas density for the AB runs, and 
between 0.01 to 100 for the BA runs. White regions indicate high density. The size of the 
box is Lx = Lz = 2rjr for run AB, and L^ = L^ = 40?7r for run BA. 
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the streaming instability can reach as high as 10'^ times background particle density. Such 
high densities are sufficient to r nake the clumps gravi t at ionally bound, promoting the forma- 



tion of planetesimals in PPDs (jJohansen et al.l 120071 ). The concentration of particles is best 



demonstrated by the cumulative probability distribution function (PDF) of particle densi- 
ties. It measures the fraction of particles whose ambient particle density exceeds a given 
value. In Figure [6] we show the particle density distribution from the series of AB and BA 
runs. Note that in the horizontal axis, we have normalized particle density to the averaged 
particle density in the simulation box. 

The ability for our code to handle the particle clump is reflected in the stiffness param- 
eter X defined in equation f lTT|) . We have tested that the maximum value of x i^i our AB 
runs is about Xmax ~ 0.2 in 64^ resolution, and it monotonically drops to Xmax ~ 0.004 in 
2048^ resolution. For the BA runs, Xmax stays less than 0.5 for most of the time and for all 
resolutions, but reaches as large as 2.0 in a few transients in the highest resolution run. The 
duration of the transients is so short that their contribution to the PDF is far below 10"^ 
and is not visible in our Figure [61 These facts show that the particle clumps are properly 
handled in our code and the obtained PDFs are not affected by the possible stiffness in the 
dense clumps. 

Our results agree quantitatively with Figure 11 of JY07 (our AB9 and BA9 256^ runs 
resolution with Np^ = 25). The PDFs calculated from our AB runs are very robust, in the 
sense that using longer period of averaging or using different initial random seeds do not 
produce any visible changes in the plot. For the BA runs, there can be some uncertainties in 
the calculation of the PDFs, mostly because that there are only a few dense clumps in the 
simulation box, the long term stochastic evolution of these big clumps makes the calculated 
PDFs more or less dependent on the period of time averaging. Nevertheless, we have averaged 
the PDFs over as long as nearly 200 orbits, which is about 170 times the correlation time of 
the clumps (which is about 7^2^^, see Figure 13 of JY07), such uncertainty is expected to be 
small and does not affect the main features to be discussed below. 

The left panels of Figure [6] show the effect of grid resolution on the density PDFs: 
higher grid resolution leads to stronger particle clumping. Both the maximum particle 
density is higher, and the number of particles residing in dense clumps is larger at higher 
resolutions. For the AB runs, the PDFs from the 512^ resolution run is almost identical 
to higher resolution runs, indicating convergence. For the BA runs, the small number of 
dense clumps in the simulation box makes convergence more difficult by averaging over a 
finite time period. Nevertheless, viewing from the PDF curves, it appears that convergence 
is finally reached at 1024^ resolution. 

We have also studied the numerical convergence with number of particles in the sim- 
ulation box. On the right panels of Figure [6|, we see that for both runs AB and BA, the 
PDFs depend very weakly on Np^. In the AB runs, the five curves almost overlap with each 
other. In BA runs, however, there is larger scatter. Again, such scatter is most likely due 
to the small number of clumps. In fact, we have compared the BA run PDFs between short 
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Fig. 6. — The cumulative probability distribution function of particle density. Top and 
bottom panels show results from the AB and BA runs respectively. On the left panels, 
number of particles per cell A'pc is fixed to 9. Green, magenta, red, blue and black curves 
label results from 64^, 128^, 256^, 512^ and 1024^ grid resolutions respectively. The dash- 
dotted curve shows the results from 2048^ resolution run for comparison. On the right panels, 
grid resolution is fixed to be 256^, while iVpc varies to be 1, 4, 9, 16 and 25, as labeled by 
green, magenta, red, blue and black curves. 
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[Te — Ts = 500r2~^) and long [Te — Tg = 1200^2"'^, shown in the figure) averaging period. The 
PDFs from longer averaging converges show better convergencqlj- We expect the scatter to 
become smaller if we run BA series for much longer time. In all, we conclude that one can 
use as small as only one particle per cell to accurately capture the density distributions of 
the solids. 



5.3. Turbulence Properties 

The streaming instability generates turbulence which differs significantly from the lam- 
inar state described by the NSH equilibrium. The turbulence is accompanied by particle 
concentration, as discussed in the previous subsection. Understanding the turbulent state 
is important since it may stro ngly affect the settli ng process as well as the radial transport 



of solids in the PPDs (JY07, iBai fc Stonell2010al ). Similar to JY07, we have computed a 



number of physical quantities to characterize the turbulence properties generated from the 
streaming instability, including: 1) The turbulent Mach number (Ma) in three directions, 
calculated from the root mean square of gas velocity fluctuations; 2) Mean radial drift ve- 
locity of the particles, normalized by the NSH radial drift velocity v^/v^'^^^', 3) Turbulent 
diffusion coefficient for particles D, calculated using the method outlined in §5.3 of JY07; 
and 4) The time and spatial averaged Reynolds stress J-" = pgUx{uy — vk), divided by the 
Reynolds stress in the NSH equilibrium (note the quantity we measure is different from J^c,x 
in JY07). Properties (1), (2) and (4) are calculated directly from spatial and time averag- 
ing from the entire simulation box. The la uncertainties are estimated from the standard 
deviation in the time sequence. In the calculation of the diffusion coefficient, we measure 
the standard deviation of the distance traveled by tracer particles at different time intervals 
ax{At),az{At). We then fit the diffusion coefficient using cr^^(At) = 2Dx,z^t. The maxi- 
mum At has been chosen to be 32^2"^ for AB runs, 320^2"^ for BA runs. The uncertainty 
of the diffusion coefficient is directly given by the linear regression estimate. The results are 
summarized in Table [H 

From Table [1] we see that the gas flow properties, namely the turbulent Mach number 
and the Reynolds stress, depend relatively weakly on the grid resolution and number of 
particles in the simulation box. This trend is also true for the radial drift velocity of the 
particles. In the AB series of runs, all these quantities agree with each other once the grid 
resolution reaches 256^. Smaller resolution runs produce somewhat different values, but with 
larger fluctuations. Note that in Figure O the typical size of particle stripes from the 1024^ 
run (AB5) is apparently smaller than that from the 256^ run (AB3). However, the statistical 
properties of gas flow from these two runs are indistinguishable. In this sense, 256^ resolution 
is sufficient to capture the essential turbulent properties for the AB runs. For BA run series. 



^° In particular, the PDF from our run BA3 (A^pc = 9) appears to have substantial lower peak densities 
than the runs with other A^pcS when we average over shorter time period. 
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the measurements from even the lowest resolution run 64^ agree very well with other higher 
resolution runs. Such numerical convergence behavior is very different from the convergence 
of particle concentration properties discussed in §5.21 

The measurements of the diffusion coefficient in the radial and vertical directions, how- 
ever, show larger variations among different runs for both AB and BA series. Such variations 
show up not only between different grid resolutions, but also at fixed resolution, varying the 
number of particles in the simulation box also makes a difference. The differences even ex- 
ceed the la uncertainties in a number of cases. Moreover, there is no systematic trend on 
such variances, especially in the BA run series. In the AB runs, it appears that towards 
higher resolution, the diffusion coefficient becomes smaller. Given the fact that bulk flow 
properties converge well at modest resolutions, the variation of diffusion coefficient measured 
from different runs, especially those with grid resolution no less than 256^, may be taken as 
uncertainties from the numerical simulations. The uncertainties in the measured diffusion 
coefficient would be about 20%, for both AB and BA runs. 



Table 1: Turbulence Properties. 



Run Resolution A^, 



pc 



Ma, 



May 1 



Ma^ ^ 



w^l^^NSH 2 



9 0.97(15) X 10-2 2.44(15) x 10' 



15) X 10-2 
07) X 10-2 
;04) X 10-2 
;03) X 10-2 

;o3) X 10-2 

05) X 10-2 
05) X 10-2 

;o4) X 10-2 
;o6) X 10-2 
;o5) X 10-2 



Vx/V:, 



D. 



D, 



77 / ttNSH 4 



ABl 
AB2 
AB3 
AB4 
AB5 
AB6 
AB7 
AB3 
AB8 
AB9 



642 

1282 

2562 

5122 

10242 

2562 

2562 

2562 

2562 

2562 



9 1.23(07) X 10-2 

9 1.24(04) X 10-2 

9 1.25(03) X 10-2 

9 1.24(03) X 10-2 

1 1.23(05) X 10-2 

4 1.24(05) X 10-2 

9 1.24(04) X 10-2 

16 1.26(06) X 10-2 

25 1.26(05) X 10-2 



2.45(07) X 10-2 
2.47(04) X 10-2 
2.47(03) X 10-2 
2.47(03) X 10-2 
2.46(05) X 10-2 
2.46(05) X 10-2 
2.47(04) X 10-2 
2.47(06) X 10-2 
2.46(05) X 10-2 



0.81 
0.98 
0.78 
0.78 
0.76 
0.76 
0.80 
0.78 
0.81 
0.83 



9 1.19(08) X 10-2 1.84(08) x 10 



1.40(80) 
1.88(27) 
2.14(07) 
2.21(06) 
2.22(05) 
2.13(08) 
2.13(08) 
2.14(07) 
2.16(09) 
2.15(08) 



2.64(08) 
4.61(05) 
5.12(04) 
4.58(04) 
4.37(06) 
5.11(06) 
5.16(03) 
5.12(04) 
5.14(04) 
5.38(05) 



X 10"^ 
X 10-^ 
X 10-5 
X 10-5 
X 10-5 
X 10-5 
X 10-5 
X 10-5 
X 10-5 
X 10-5 



7.01(29) 
9.50(63) 
3.01(06) 
2.89(07) 
2.28(02) 
2.76(07) 
3.15(10) 
3.01(06) 
3.29(07) 
3.25(08) 



X 10-5 
X 10-5 
X 10-5 
X 10-5 
X 10-5 
X 10-5 

X 10-5 
X 10-5 
X 10-5 
X 10-5 



X 10-3 
X 10-3 
X 10-3 

X 10-3 
X 10-3 
X 10-3 
X 10-3 
X 10-3 
X 10-3 
X 10-3 



1.72(79 
2.38(30 
2.67(12 
2.69(10 
2.68(09 
2.67(16 
2.67(15 
2.67(12 
2.71(16 
2.72(16 



BAl 
BA2 
BA3 
BA4 
BA5 
BA6 
BA7 
BA3 
BA8 
BA9 



642 

1282 

2562 

5122 

10242 

2562 

2562 

2562 

2562 

2562 



9 

9 

9 

9 

1 

4 

9 

16 

25 



1.20(10) X 10-2 
1.21(11) X 10-2 
1.16(11) X 10-2 
1.36(09) X 10-2 
1.17(12) X 10-2 
1.12(05) X 10-2 
1.21(11) X 10-2 
1.16(10) X 10-2 
1.14(08) X 10-2 



-2 

1.87(10) X 10-2 
1.81(11) X 10-2 
1.83(11) X 10-2 
1.88(09) X 10-2 
1.86(12) X 10-2 
1.80(05) X 10-2 
1.81(11) X 10-2 
1.83(10) X 10-2 
1.81(08) X 10-2 



3.88 
4.03 
3.88 
3.93 
3.89 
4.00 
3.92 
3.88 
4.00 
3.89 



08) X 10-2 

10) X 10-2 

11) X 10-2 

11) X 10-2 
;09) X 10-2 

12) X 10-2 
05) X 10-2 

L) X 10-2 
10) X 10-2 
;08) X 10-2 



0.65(07) 
0.65(06) 
0.66(07) 
0.70(06) 
0.60(06) 
0.67(08) 
0.71(04) 
0.66(07) 
0.69(05) 
0.70(06) 



2.18(06) 
2.36(06) 
2.30(09) 
2.29(11) 
2.20(05) 
2.71(11) 
1.87(03) 
2.30(09) 
2.15(06) 
2.23(08) 



1.01(04) 
1.37(06) 
0.97(05) 
1.39(05) 
1.16(08) 
1.44(05) 
1.23(05) 
0.97(05) 
1.35(06) 
1.32(05) 



X 10-2 
X 10-2 
X 10-2 
X 10-2 
X 10-2 
X 10-2 

X 10-2 
X 10-2 
X 10-2 
X 10-2 



0.5^7 
0.61'(06 
0.62(08 
0.67(06 
0.56(06 
0.64(08 
0.68(04 
0.62(08 
0.66(05 
0.67(06 



The number in parenthesis quotes the la uncertainty of the last two digits. See §5.31 for details. 
^ Mach number in the radial, azimuthal and vertical directions. 

2 Radial drift velocity of particles, normalized by the NSH value. 

3 Turbulent diffusion coefficient of the particles in the radial and vertical directions. 
^ Mean Reynolds stress of the gas, normalized by the NSH value. 
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6. Summary and Conclusion 

We have presented the implementation of a hybrid particle-gas scheme in the grid-based 
Athena MHD code. The particle and gas are assumed to be coupled aerodynamically, and 
a size distribution of particle species with different stopping times is allowed. Our imple- 
mentation is extendable to include gravitational coupling, which is left for future work. The 
main purpose for the code development is to study the dynamics of gas and solids in proto- 
planetary disks (PPDs). The solid size range where the aerodynamic coupling has significant 
effects in PPDs roughly spans from millimeter to a few tens or hundreds meter sized bodies, 
depending on the disk mod el and location, and this is the regime most relevant for study- 



ing planetesimal formation (jChiang fc Youdinll2009l ). In this pape r we mainly address the 



numerical method and code performance. In a forthcoming paper (JBai fc Stondl2010al ). we 
will describe applications to PPDs. 

The numerical algorithm of our hybrid code is based on a second-order accu rate predictor- 



corrector scheme. The algorithm is different from other existing codes (YJOT jBalsara et al. 



20091 : lMiniatill2010l ). and is very simple and robust. Our implementation of particle-gas cou- 
pling is fully conservative: backreaction from the particles to the gas is treated carefully that 
the total momentum is conserved exactly. Our hybrid code works well in non-stiff regime of 
particle-gas coupling, as well as the stiff regime without significant particle mass loading [i.e., 
the parameter x defined in equation (ITT]) does not exceed order unity]. This is made possible 
by two implicit particle integrators. These include a semi-implicit integrator, generalized 
from the second order leap-frog integrator, which has much better stability properties than 
any explicit methods, and which preserves the geometric properties of particle orbit exactly 
in the limit of zero drag force. In addition, a fully-implicit integrator is designed for treating 
extremely stiff problem when particle stopping time is much smaller than the simulation 
time step. We have extensively tested our code performance, including the linear growth 
rate test of the streaming instability. Subsequent test on the non-linear saturation of the 
streaming instability further confirms that our code performs as good as the higher order 
finite difference Pencil code, as the flow properties measured from our simulations agree well 
with the results by JY07. 

We have also studied the numerical convergence of our method in the non-linear regime 
of the streaming instability. We pick two representative runs from JY07 and have performed 
a series of simulations by varying grid resolutions and total number of particles in the simu- 
lation box. We find that the convergence properties in the non-linear regime is very different 
and more complicated than those in the linear regime. The main conclusions drawn from 
this convergence study are summarized below. 

1. Requirement for numerical convergence strongly depends on particle stopping time r^. 
For Ts = 0.1, about 128 cells per rjr is necessary for convergence. 

2. Equal number of particles and grid cells is sufficient for numerical convergence. 
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3. Gas flow properties converge very well at modest grid resolution and is not sensitive 
to number of particles used in the simulation. 

4. Particle concentration properties converge at very high grid resolution. Higher resolu- 
tion leads to stronger clumping. 

5. Particle diffusion properties depend on the numerical setup in a subtle way, leading to 
about ±20% uncertainties in the measurements. 

These convergence tests provide useful information on the reliability and uncertainty of 
this kind of hybrid simulations of the particle-gas interaction. They also provide a guide on 
the choice of grid resolution and number of particles for future simulations of similar and 
more complicated problems. Although all the results are based on vertically unstratified 
simulations, we may generalize these criteria to simulations with vertical gravity. Also, with 
more than one particle species, we expect numerical convergence with one particle per cell 
per particle species, where different particle species have different stopping times. 



Recently, iRein et al.l (120 lOf ) addressed the validity of the super-particle approximation 
in numerical simulations with particles. They emphasized that it is essential to maintain 
the important timescales in the scaled system (i.e. the simulation) to be equivalent to the 
timescales in the real system, otherwise one cannot achieve numerical convergence. This 
is unlikely to be relevant to the problem considered in this paper because we have not 
added gravitational interaction nor particle collisions. The only time scales to be fixed 
are the orbital period and the stopping time, which are combined into the dimensionless 
stopping time Tg. Our results imply that numerical convergence can still be achieved when 
aerodynamic interaction between particle and gas is added, strengthen conclusions of Rein 
et al. 

We thank Anders Johansen and Andrew Youdin for helpful discussions in the "Dynamics 
of Discs and Planets" programme at the Newton Institute, Cambridge University. We are 
grateful to Paul Bode and Anatoly Spitkovsky for advises on the data structure of the particle 
component, and Emmanuel Jacquet for his earlier development of the particle module. We 
thank our referee for carefully reading the manuscript and the helpful suggestions. This work 
is supported by NSF grant AST-0908269. XNB acknowledges support from NASA Earth 
and Space Science Fellowship. 



A. Additional Linear Growth Rate Test of the Streaming Instability 

The linear and non-linear SI simulations in this paper adopt relatively large particle 
stopping times r^ < 0.1. In this appendix we provide two additional linear growth test 
problems to test code performance for smaller r^. The simulation set up is described section 
14. 3i the first test ("linC"), r^ = 0.01, e = 2, K^ = K^ = 1500, the analytical growth rate is 
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s = 0.5981. In the second test ("linD"), r^ = 10-^ e = 2, and K^ = K;, = 2000, with growth 
rate s = 0.3154. Both modes are chosen to be close to the fastest growth mode. However, 
these tests are more computationally costly because of the larger K. The eigenvectors of the 
two modes are listed in Table Htj. The table format is the same as Table 1 of YJ07. 

In Figure [7] we show the numerical growth rate measured as a function of grid resolution 
for the two tests "linC" and "linD". The CFL number is fixed at 0.4 for these tests. Note 
that for smaller r^, the perturbations on the gas is extremely small, and capturing the 
correct growth rate on the gas is very difficult. For other quantities, numerical convergence 
to theoretical growth rate is clearly reached as shown in the figure, however, for smaller r^, 
higher resolution of about 128 cells per wavelength is needed for numerical convergence. 



'^We thank Andrew Youdin for providing the eigenvectors. 
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Fig. 7. — Measured growth rate of the streaming instabihty eigenmodes "hnC" (diamonds) 
and "hnD" (asterisks) using the semi-imphcit integrator. Sohd hne marks the theoretical 
growth rate, and dashed hues with symbols label the measured growth rate as a function of 
grid resolution (number of cells per wavelength). We use CFL= 0.4 in all the tests. 



Table 2: Test Mode Eigensystems. 



Test 



M.x 



M„ 



U^, 



Pa 



Wr 



W,. 



w. 



OJ 



linC Ts = 10^2^ e = 2 
{K^ = K, = 1500) 

linD r, = 10-^, e = 2 
(/r, = K, = 2000) 



-0.1598751 0.1164423 

-0.0079669i +0.0122377i 
-0.1719650 0.1918893 

-0.0740712i +0.0786519i 



0.1598751 8.684872e-8 

-0.0079669i +5.350037e-7i 

0.1719650 2.954631e-7 

-0.0740712i +1.141385e-7i 



-0.1567174 0.1159782 0.1590095 0.1048236 

-0.0028837i +0.0161145i -0.0024850i +0.5980690i 

-0.1715840 0.1918542 0.1719675 0.3224884 

-0.0740738i +0.078737H -0.0739160i +0.3154373i 



Note: the table format is the same as Table 1 of YJ07. 
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